//QUIRKSIM
//Copyright (c) 2008 Gregory Andrew Robinson (garobinson@ucdavis.edu)
//
//Permission is hereby granted, free of charge, to any person
//obtaining a copy of this software and associated documentation
//files (the "Software"), to deal in the Software without
//restriction, including without limitation the rights to use,
//copy, modify, merge, publish, distribute, sublicense, and/or sell
//copies of the Software, and to permit persons to whom the
//Software is furnished to do so, subject to the following
//conditions:
//
//The above copyright notice and this permission notice shall be
//included in all copies or substantial portions of the Software.
//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
//
//EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
//OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
//NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
//HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
//WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
//OTHER DEALINGS IN THE SOFTWARE.

#include <iostream>
#include <fstream>
#include <gflags/gflags.h>
#include <glog/logging.h>
#include "particle.h"

using quirksim::particle::QuirkPair;
using quirksim::particle::MassiveParticle;

DEFINE_double(timestep, 1.0e-2,
              "Time interval between numerical integration iterations.");

namespace particle {
  double c = 1.0;
}

int main(int argc, char* argv[]) {
  double time;
  long percent;
  long i;
  long num_iter = 10000;
  double mass = 10e-28;
  double sigma = 1000;

  Hep3Vector relative_p;

  google::ParseCommandLineFlags(&argc, &argv, true);
  google::InitGoogleLogging(argv[0]);

  std::cout.precision(20);

  Hep3Vector position1(-1, 0, 0);
  Hep3Vector momentum1(-1.2239283748106613, 1.2239283748106613, 0);
  Hep3Vector position2(1, 0, 0);
  Hep3Vector momentum2(1.2239283748106613, 1.2239283748106613, 0);

  MassiveParticle quirk1(mass, position1, momentum1);
  MassiveParticle quirk2(mass, position2, momentum2);
  quirksim::particle::QuirkPair pair(quirk1, quirk2, sigma);

  relative_p = quirk1.position() - quirk2.position();

  std::ofstream fout("quirkoutput.csv");

  if(!fout.is_open()) {
    std::cout << "could not open file.";
    return 1;
  }

  for(i=0; i<=num_iter; i++) {
      time = i*FLAGS_timestep;
      percent = i*100L/num_iter;
      LOG(INFO) << "t=" << time << " seconds (" << percent
                << "%):\n";
      LOG(INFO) << "q1 x: "
                << quirk1.position().x() << ", "
                << quirk1.position().y() << ", "
                << quirk1.position().z();
      LOG(INFO) << "q2 x: "
                << quirk2.position().x() << ", "
                << quirk2.position().y() << ", "
                << quirk2.position().z() << "\n";
      LOG(INFO) << "q1 p: "
                << quirk1.momentum().x() << ", "
                << quirk1.momentum().y() << ", "
                << quirk1.momentum().z();
      LOG(INFO) << "q2 p: "
                << quirk2.momentum().x() << ", "
                << quirk2.momentum().y() << ", "
                << quirk2.momentum().z() << "\n";
      fout << time << ", "
           << quirk1.position().x() << ", "
           << quirk1.position().y() << ", "
           << quirk1.position().z() << ", "
           << quirk2.position().x() << ", "
           << quirk2.position().y() << ", "
           << quirk2.position().z() << ", "
           << quirk1.momentum().x() << ", "
           << quirk1.momentum().y() << ", "
           << quirk1.momentum().z() << ", "
           << quirk2.momentum().x() << ", "
           << quirk2.momentum().y() << ", "
           << quirk2.momentum().z() << "\n";
    try{
      pair.Evolve(FLAGS_timestep);
    } catch (std::exception e) {
      std::cout << e.what() << " ...wtf\n";
    }
  }
  fout << std::flush;
  fout.close();
  return 0;
}
